%PCA on FlowJo raw data for clustered gate
%12/11/12 Thomas Stevens
%Input: Cluster Q2 FITC, B as CSV
%% FOR EPS, need to use painters renderer (shading from OPENGL destroyed, 
%% can add back in Illustrator) and CMYK colors

data_equal = csvread('../All_Lim_FACS_tvs/PCA_Raw_FITC_B_111112_t6h30m/1-1_maxdox.csv');
data_redbias = csvread('../All_Lim_FACS_tvs/PCA_Raw_FITC_B_111112_t6h30m/g1-r8_maxdox.csv');
%Error: Exported all columns, FITC, B are 3,4
data_greenbias = csvread('../All_Lim_FACS_tvs/PCA_Raw_FITC_B_111112_t6h30m/g8-r1_maxdox.csv');
data_greenbias(:,[1 2 5]) = [];

%log before linear components visible
%CheckPlot
dscatter(log10(data_equal(:,1)),log10(data_equal(:,2)))
figure
dscatter(log10(data_redbias(:,1)),log10(data_redbias(:,2)))
figure
dscatter(log10(data_greenbias(:,1)),log10(data_greenbias(:,2)))

%PCA - coeff is orthonormal basis for PC space
[coeff,score,latent] = princomp(log10(data_equal));
%First PC explains 89.4% of variability
cumsum(latent)./sum(latent)
%Check with biplot
figure
biplot(coeff(:,1:2),'Scores',score(:,1:2))
%Project other data onto 1:1 pc ** Need to recenter by geo mean of 1:1
center_redbias = log10(data_redbias) - ones(size(data_redbias,1),1)*log10(mean(data_equal));
center_greenbias = log10(data_greenbias) - ones(size(data_greenbias,1),1)*log10(mean(data_equal));
score_redbias = center_redbias*coeff;
score_greenbias = center_greenbias*coeff;
%Check scatter of renormalized data
figure,hold on
scatter(score_redbias(:,1),score_redbias(:,2),'r')
scatter(score_greenbias(:,1),score_greenbias(:,2),'g')
scatter(score(:,1),score(:,2),'k')

%Histograms along each pc
figure, hold on
%PC1 - size
[n1,x1] = hist(score(:,1),sqrt(numel(score(:,1))));
[n1r,x1r] = hist(score_redbias(:,1),sqrt(numel(score_redbias(:,1))));
[n1g,x1g] = hist(score_greenbias(:,1),sqrt(numel(score_greenbias(:,1))));

%PC2 - composition
figure, hold on
[n2,x2] = hist(score(:,2),sqrt(numel(score(:,2))));
[n2r,x2r] = hist(score_redbias(:,2),sqrt(numel(score_redbias(:,2))));
[n2g,x2g] = hist(score_greenbias(:,2),sqrt(numel(score_greenbias(:,2))));

%Check plots - size
figure, hold on
plot(x1,n1,'k')
plot(x1r,n1r,'r')
plot(x1g,n1g,'g')
%Check plots - composition
figure, hold on
plot(x2,n2,'k')
plot(x2r,n2r,'r')
plot(x2g,n2g,'g')

%Kernel density
N = {score score_redbias score_greenbias};
F = [];
X = [];
for j=1:2
    for i = 1:size(N,2)
        [f,xi] = ksdensity(N{i}(:,j));
        F = [F f'];
        X = [X xi'];
    end
end
%Plot trans-shaded kernel densities for each PC
figure, hold on
colors = 'brg';
for i = 1:3
    h = fill(X(:,i),F(:,i),colors(mod(i-1,3)+1));
    set(h, 'EdgeColor','none', 'FaceAlpha', 0.5);
end
defaultPlot()
xlabel('PC1 ~ Cluster Size')
ylabel('Probability')
figure, hold on
for i = 4:6
    h = fill(X(:,i),F(:,i),colors(mod(i-1,3)+1));
    set(h, 'EdgeColor','none', 'FaceAlpha', 0.5);
end
defaultPlot()
xlabel('PC2 ~ Cluster Composition')
ylabel('Probability')